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ABSTRACT 



The numerical experiments, carried out through the use of the vorticity- 
stream function equations and their finite difference form, on sinusoidally- 
oscillating as well as co-existing flows (sinusoidal oscillation plus steady mean 
flow) at low and intermediate Keulegan-Carpenter numbers are described. A 
third-order in time, second-order in space, three-level predictor-corrector finite- 
difference scheme has been used. The Poisson equation for the stream function 
was solved by a Fast Poisson Solver based on the High Order Difference 
Approximation with Identity Expansion (HODIE) and the Fast Fourier 
Transform (FFT) methods provided by the National Center for Atmospheric 
Research for the solution of separable elliptic partial differential equations with a 
non-square grid. The analysis has produced force-transfer and fluid-damping 
coefficients comparable to those obtained experimentally for both types of flows 
(i.e., with and without current) and to those obtained with a square grid through 
the use of the IMSL library. 
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NOMENCLATURE 



A = amplitude of sinusoidal oscillations 

a = transformation parameter 

Cjl = inline force coefficient 

Cl = transverse force coefficient 

Cjj = Fourier-averaged drag coefficient 

C m = inertia coefficient = 1 + k'jj 

D = diameter 

Fil = inline force 

Fp = transverse force 

K = Keulegan-Carpenter Number = U m T/D 

K+ = K by current = K(1 + U 0 /U m ) = K + U 0 T/D 

= Fourier-averaged added mass coefficient 

n = order of polynomial (see Eq 17) 

p s = pressure on cylinder 

Poo = pressure at outer boundary 

p s = pressure coefficient = (p s - p^) / pUoo 
R = cylinder radius 

Re = Reynolds number VU/v, or UmD/v 



IX 



Re + = Re modified by current = Re(l + Uo/Um) = Re + UoD/v 
r = radial distance 

T = period of oscillations 

t = time 

U = time dependent velocity 

Um = maximum velocity in pure sinusoidal flow 

Uo = collinear steady current velocity 

V = constant velocity reached at the end of the acceleration period 

V r = current ratio = Uo/Um 

(3 = D 2 /vT = Re/K 

Aa = disturbance oscillation applied to the ambient flow, (in degrees) 

AS, = computational grid spacing 

Q = relative displacement of fluid during the acceleration period 
= (S/R) v = 0.5 (dU/dt) t l/R = 0.5 V t v /(R) 

|! = dynamic viscosity 

v = kinematic viscosity 

p = density 

0 = angular position 

\f / = stream function 

(0 = vorticity 
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I. INTRODUCTION 



Numerical experimentation in fluid dynamics, through the use of 
finite-difference, finite-element, and discrete-vortex methods, has attracted 
considerable attention during the past two decades and produced laminar 
flows difficult to measure and turbulent flows hard to verify and impossible 
to generalize. The reasons for this are relatively simple. Numerical solutions 
based on the full Navier-Stokes equations are not stable at high Reynolds 
numbers and the instability is non-linearly related to the particular flow, 
input parameters and the discretization conditions. Also, the real flow at the 
computed Reynolds numbers may be turbulent, at least in some regions of 
the flow, and the numerical experiment does not imitate the physical 
experiment. Furthermore, the observed physical and numerical instabilities 
do not necessarily correspond to each other. Assuming that the calculations 
for a given flow are carried out at sufficiently small Reynolds numbers, where 
the flow is known to remain stable and laminar, one quickly discovers that it 
is practically impossible to measure, to any credible degree of accuracy, most 
or all of the predicted quantities (except the Strouhal number and the 
photographs of the flow patterns). 

Evidently, one’s view of the state of the numerical modelling depends to 
a large extent on one's objectives. For example, if the objective is to obtain 
some approximate answers and flow kinematics, one might be perfectly 
satisfied with the existing codes. If the objective is to match the measured and 
calculated results (e.g., lift and drag coefficients), one might achieve the 
desired objective by fine tuning a number of model parameters (e.g., the order 
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of approximation of the velocity and/or vorticity gradients, particularly near 
the wall, mesh size, time step, type of discretization, outer boundary, just to 
name a few). If one's objectives are to perform numerical experiments for 
sake of numerical experiments, with no concern with the compatibility of the 
numerical and experimental results, then one can objectively asses the model 
instead of attempting to attribute to it artificial powers of prediction. 

As far as the turbulent flows are concerned, some or all of the predictions 
of the numerical calculations for a given flow depend on the closure model 
used. Some models do better than others for some flows and worse than 
others for other flows. No model, however sophisticated, has a corner on the 
numerical market. Among the numerous theoretical, numerical, and 
experimental investigations, impulsively-started steady flow about a circular 
cylinder has occupied a prominent place partly because of its intrinsic interest 
towards the understanding of the evolution of separation, vortex formation, 
growth, and partly because it provided the most fundamental case for the 
comparison and validation of various numerical methods and codes. In 
recent years, attention has turned to a broader class of relatively manageable 
time-dependent flows about bluff bodies: Non-impulsively-started flows, 
sinusoidally oscillating flows, co-existing flows (uniform flow plus oscillating 
flow), flow from one steady state to another (at a lower or higher Reynolds 
number through the use of prescribed changes in velocity), and so on. The 
solution of these problems at sufficiently high Reynolds numbers will have 
far reaching theoretical and practical consequences. As noted above, this is 
not yet the case, and the solutions must necessarily be confined to cases where 
the accurate prediction of physical experiments and the instant gratification 
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are not the real objectives. However, it is hoped that even the approximate 
solutions will have enough information to elucidate the physics of the 
phenomenon. 

The hydrodynamic loading situations which are well understood are 
those which do not involve flow separation. Thus, they are amenable to 
nearly exact analytical treatment. These concern primarily the determination 
of the fluid forces on large objects in the diffraction regime where the 
characteristic dimension of the body relative to the wave length is larger than 
about 0.2. The use of various numerical techniques is sufficient to predict 
accurately the forces and moments acting on the body, provided that the 
viscous effects and the effects of separation for bodies with sharp edges are 
ignored as secondary. 

The understanding of the fluid-structure interactions which involve 
extensive flow separation and dependence on numerous parameters such as 
Reynolds number, Keulegan-Carpenter number, relative roughness, relative 
motion of the body, proximity effects, hydroelastic response, etc. is far from 
complete (Sarpkaya & Isaacson 1981). There are several reasons for this. First, 
although the physical laws governing the motion (the Navier-Stokes 
equations) are well understood, valid approximations necessary for 
numerical and physical model studies are still unknown. Even the 
unidirectional steady flow about a bluff body remains theoretically 
unresolved. Much of our understanding of vortex shedding behind bluff 
bodies came from steady-flow experiments, highly idealized models, and 
limited numerical solutions. Most of the numerical studies based on the use 
of the Navier-Stokes equations and some suitable spatial and temporal 
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differencing schemes are limited, out of necessity, to low Reynolds number 
flows. A second reason why progress has been slow is that the bluff body 
problems involving wake return are an order of magnitude more complex 
and there has been only a handful of limited applications of the methods 
based on Navier-Stokes equations. 

The formation of a wake gives rise not only to a form drag, as it would be 
the case if the motion were steady, but also to significant changes in the 
inertial forces. The velocity-dependent form drag is not the same as that for 
the steady flow of a viscous fluid, and the acceleration-dependent inertial 
resistance is not the same as that for an unseparated unsteady flow of an 
inviscid fluid. In other words, the drag and inertial forces are interdependent 
as well as time-dependent. These effects are further compounded by the 
diffusion and decay of vortices and by the three-dimensional nature of 
vorticity due to turbulent mixing, finite spanwise coherence, and the random 
nature of the vortices (which give rise to cycle-to-cycle variations and 
numerous flow modes even under controlled laboratory conditions). The 
stronger and better correlated the returning vortices, the sharper and more 
pronounced the changes are in pressure distribution on the body and in the 
integrated quantities such as the lift, drag, and inertia coefficients. 

It is clear from the foregoing that the objectives of the present 
investigation are to carry out extensive numerical experiments through the 
use of the vorticity-stream function form of the Navier-Stokes equations and 
their finite difference form, on co-existing flows (sinusoidal oscillation plus 
steady mean flow). The expectations are that the results will point out the 
strengths and weaknesses of the code, for the particular type of formulation 
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used, explain the reasons between the various numerical predictions of the 
same problem, and, hopefully, shed some light on the physics of flows 
heretofore uncalculated. 
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II. BACKGROUND STUDIES 



A finite difference analysis of the Navier-Stokes equations for a 
sinusoidally-oscillating ambient flow about a circular cylinder at K (Keulegan- 
Carpenter Number) = U m T/D = 5 (Re = 1000) and K = 7 (Re = 700) has been 
attempted by Baba & Miyata (1987). Their results have shown that the 
calculations can be carried out only for short times (less than two cycles of 
flow oscillation) with a non-super computer. Murashige, Hinatsu and 
Kinoshita (1989) have used a similar method to analyze three cases (K = 5, 7, 
and 10) at higher Reynolds numbers around 10 4 . The flow was perturbed by 
artificial means to trigger an asymmetry. At K = 10, a transverse vortex street 
appeared, in agreement with experimental observations. The numerical 
simulation of steady flow past a circular cylinder undergoing in-line and/or 
transverse oscillations through the use of two-dimensional unsteady Navier- 
Stokes equations was undertaken by Lecointe et al. (1987) for relatively small 
amplitudes (A/D = 0.13). Justesen (1991) presented extensive results obtained 
from a numerical solution of a vorticity-stream function formulation of the 
Navier-Stokes equations for the flow around a circular cylinder in planar 
oscillating flow at small Keulegan-Carpenter numbers in the subcritical 
Reynolds number range. Justesen introduced a straining parameter "a" in 
order to better resolve the large gradients near the cylinder surface. This is in 
addition to the logarithmic straining, commonly used as part of the 
transformations, for a better resolution of the gradients near the body. 
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Evidently, Justesen's transformation for a = 0 defaults to the logarithmic 
straining. However, "a" becomes another disposable parameter, dependent 
on at least K and Re. Justesen had to choose judiciously the value of the 
straining parameter for each K in order achieve drag and inertia coefficients 
in satisfactory agreement with those obtained experimentally. A systematic 
numerical variation of the governing parameters for an arbitrary U(t) is 
extremely difficult. 

The in-line oscillations of a cylinder in uniform flow (or the sinusoidally 
oscillating flow with a steady mean flow) has been the subject of intense 
interest in recent years (see, e.g., Sarpkaya & Isaacson, 1981 and Sarpkaya & 
Storm, 1985) in connection with the understanding of the behavior of hot- 
wire anemometers and the fluid loading of structures subjected to currents, 
gusts and other types of unsteady flows. The biassing of the shedding of the 
vortices by the current causes profound changes in both the drag and inertia 
coefficients, relative to their no-current values. The mobile separation points 
undergo large excursions, as much as 120 degrees during a given cycle of 
oscillation over a circular cylinder (Sarpkaya and Butterworth, 1992). These 
effects are further compounded by the diffusion and decay of vortices and by 
the three-dimensional nature of vorticity due to turbulent mixing, reduced 
spanwise coherence, mutual-induction instability, and the random nature of 
vortices which give rise to cycle-to-cycle variations and numerous flow 
modes even under controlled laboratory conditions. It is because of these 
reasons that the present work is dedicated to the numerical analysis of such 
flows. 
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III. NUMERICAL REPRESENTATION 



A. COMPUTATIONAL METHOD 

Here only a brief description of the computational method is presented. 
A more in depth description is given by Wang (1989) and Putzig (1991). 

The fluid is assumed to be two-dimensional, incompressible and viscous. 
The governing equations for the solution are the Navier-Stokes equations, 
with the stream function and the vorticity as independent variables. To 
achieve a higher density of mesh points near the cylinder surface, the 
computational domain is transformed from the physical plane in Fig. 1 to a 
rectangular plane in Fig. 2 (see Appendix). In the rectangular plane, the mesh 
is maintained at a uniform grid spacing. It is necessary to have more mesh 
points closer to the cylinder surface because in this region the gradients of 
both the vorticity and the stream function are the largest. 

A third-order in time, second-order in space, three-level predictor- 
corrector finite-difference scheme is used to solve the vorticity-transport 
equation. A Fast Poisson Solver automatically discretizes the separable 
elliptic equation which is then solved by a generalized cyclic reduction 
algorithm. One of the solvers was provided by the IMSL mathematics library, 
for the solution of separable elliptic partial differential equations with a 
square grid, and the other by the National Center for Atmospheric Research, 
for the solution of separable elliptic partial differential equations with a non- 
square grid. 
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The unsteady Navier-Stokes equations in the polar coordinates, as 
defined by the vorticity transport equation and the vorticity/stream-function 
equation are, 



and 



00) 1 
0t r 



dr ae ae ar 



= vV 2 co 



(1) 



where 



V 2 y = o 



^2 d 2 i a i a 2 

= — =•+— 



(2) 



0r 2 r0r a0 ^ 



(3) 



co and \\f are the vorticity and the stream function, v is the kinematic viscosity, 
t is the time and, r and 0 are polar coordinates directions (see Fig. 1). The 
velocity components in the r and 0 directions are defined by 



1 0W 011/ 

u = — v = — L 

r 09 and 0r 

The boundary conditions for the physical problem are: 

(1) no slip and zero normal velocity on the surface of the cylinder 

\| f = — = 0 on r = R 

0r 

and (2) the potential flow at infinity is defined as 



(4) 



(5) 



V = U(r--)sin0 
r 



(6) 



and to = 0 at r = U is the external flow and R is the radius of the cylinder. 
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The coordinate transformations required to go from the physical domain 
to the computational domain are: 



The transformation of the non-dimensionalized vorticity-stream 
function equations and their finite difference form through the use of the 
central difference approximation for vorticity and a two-step, three-level, 
predictor-corrector scheme, with a third order accuracy in time, are described 
in detail in Wang (1989), in Fredrickson (1990) and in Putzig (1991) and will 
not be repeated here. 

B. CALCULATION OF THE FORCE COEFFICIENTS 

The in-line and transverse force coefficients are determined from the 
combined contributions of the shear and pressure forces acting on the 
cylinder. The viscous forces are calculated from t s = |i(D. The total in-line 
force then reduces to 



r = R exp(a^) and 9 = aq 



(7) 



where R is the radius of the cylinder and 'a' is a transformation parameter. 



F 1L = — jp 11 p s cos(6)Rde - Jg " (0.0) sin(w)Rd6 



( 8 ) 



and the total lift force as 



Fl = -J 0 27t p s sin(6)Rd0 - jQ 7I |icocos(co)Rd0 



(9) 



After dividing the in-line and the lift-force equations by (0.5 pU 2 D) and defining 
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( 10 ) 



= _ 2<Ps - P~) 
Ps - tt 2 

pU^ 



the force coefficients reduce to 



c il = --Jo"p s cos( 8 )Rd 6 - — } 0 2 "pcosin(co)Rd 8 
2 Re 



( 11 ) 



and 



c l = -1 Jo” p s sin( 6 )Rd 9 - ^ J o 2 ,I pcocos(co)Rd 0 



(12) 



The pressure coefficient is determined from the Navier-Stokes 
equations in terms of dimensionless vorticity. Once integrated with respect to 
0, one has 



Ps( 0 ) = Ps(°) + — t 



4 r2rcf0aO 



Re J 0 ^ 5r y r=1 



d0 



( 13 ) 



Equation (13) is substituted into equations (11) and (12) to determine the 
numerical scheme for the in-line and transverse force coefficients. 



q = — r 

^ i?Jo 



H 




cos(0) + o3sin(0) ► 


Jo 


1 

\ 

o 

1 


j 



d0 



( 14 ) 



and 



C L = 



Re Jo 



rj 


— i 

rS 

i SI 


d0 


sin(0) + O3Cos(0) ► 


Jo 


dr ) 


II 

>— i 





de 



( 15 ) 



11 



The radial derivative of the vorticity on the surface of the cylinder, 
appearing in Equations (14) and (15), is determined through the use of discrete 
pointwise approximations of various orders, ranging from second to tenth 
order. For a second order approximation, one has 



(M) _ -3Q3i + 4(O i+ - 1 -03 i+2 f Q( A ^2) 

dr J r=1 2A£, 



(16) 



For higher order polynomials Equation (16) may be written as 



( 5(0i ^ _ AcOj + Bco i+1 + C(o i+2 + Dco i+3 + Eco i+4 + . 

dr J r=1 AS, 



+ 0(A^ n ) (17) 



in which the coefficients A -K are given in Table 1 below. 
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Table 1: Coefficients of the Polynomial in Eq. (17) 





n=2 


n=4 


n=6 


n=8 


n=10 


A 


-3/2 


-25/12 


-49/20 


-761/280 


-7381/2520 


B 


2 


4 


6 


8 


10 


C 


-1/2 


-3 


- 15/2 


-14 


-45/2 


D 




4/3 


20/3 


56/3 


40 


E 




-1/4 


-15/4 


-35/2 


-105/2 


F 






6/5 


56/5 


252/5 


G 






-1/6 


-14/3 


-35 


H 








8/7 


120/7 


I 








-1/8 


-45/8 


J 










10/9 


K 










-1/10 
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C CALCULATION OF THE DRAG AND INERTIA COEFFICIENTS 

If one were to associate the total force with a velocity-square-dependent 
drag force and an acceleration-dependent inertial force then the coefficient 
associated with the latter may be interpreted as some measure of the added 
mass. But one must bare in mind that such a decomposition is far from being 
unique. 

It has been customary to express the fluid force acting on a body moving 
in a fluid otherwise at rest as 

F(t) = ipC^ Ap|{(U 0 + U(t)}|{(U 0 + U(t)} + pk,H (18) 

where U G represents the steady velocity; U(t), the time-dependent 
oscillations; C^ / the Fourier-averaged drag coefficient and k~ , the Fourier- 

averaged added-mass coefficient. It is customary to use an inertia coefficient 
C m for a fluid in motion about a body at rest through the use of C m = 1 + kVj . 

The Fourier averages of the drag and added-mass coefficients over a 
period of T may be calculated by multiplying both sides of Equation (8) once 
with U(t) and once with dU/dt to yield 

T 

2jF(t)U(t)dt 

C d u = j 2 (19) 

pA p J|{(U D + U(t)}| {(U 0 + U(t)} U(t)dt 
0 
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and 



(Cm - 1) = 

pVj[dU(t)/dt] 2 dt 

0 



jF( t )^dt 
0 



(20) 



which may be evaluated readily provided that sufficiently reliable data are 
available for F(t), U G , U(t), and dU(t)/dt. 



A simple dimensional analysis of the flow under consideration 



a relative amplitude or Keulegan-Carpenter number, Mach number, 
Reynolds number, and a parameter involving Uo (e.g., UoT/D or 
Uo/[U(t)]max)- There are numerous possibilities regarding the definitions of 
the relative amplitude or Keulegan-Carpenter number and the Reynolds 
number. The purpose of the search for a more suitable Keulegan-Carpenter 
number and/or Reynolds number is to enhance the correlation of the data to 
reduce the number of the governing parameters, possibly eliminating UoT/D 
as an independent parameter. The list of possible Reynolds numbers and 
Keulegan-Carpenter numbers is long and will not be given here. Suffice it to 
note that the two force-coefficients for the flow about a cylinder may be 
written as 



shows that the time-averaged force coefficients (C^j and k^- ) are functions of 




(21) 
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or as 




= fi(K + ,Re + ,VK) 



(22) 




in which 



K = U m T/D , Re = U m D/v , VK = U 0 T/D 



K+ = K(1 + U(/U m ) = K + U 0 T/D 



Re + = Re(l + Uo/Uni) — Re + UoD/v 



(23) 



The purpose of the present calculations was not to provide a detailed 
comparison between the measured and calculated forces but rather to attempt 
to establish a relationship between the shedding of vortices and the relative 
magnitude of the current. The particular values of K and V r chosen for the 
calculations (K = 4 - 6, V r = 0.0 - 1.2) was one for which some experimental 
data were available at comparable (3, Re, and V r values. 
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IV. DISCUSSION OF RESULTS 



The numerical experiments were carried out through the use of a VAX- 
2000, a VAX-3520, a CRAY Supercomputer. The solution procedure and 
technique have been validated for several types of unsteady flows, i.e., 
impulsively-started, suddenly-stopped, and uniformly-decelerated flows 
before applying it to co-existing flows (oscillatory flow or sinusoidally- 
oscillating flow with mean velocity). Excellent agreement with flow 
visualization and experimentally determined drag and lift coefficients has 
been obtained for both symmetric and asymmetric wake solutions. It is this 
validation that led to the exploration of the characteristics of sinusoidally- 
oscillating flows superimposed on a mean velocity. The oscillation was 
specified by U = U 0 +U m sin(27it/T) in which U G is the steady mean velocity 
and U m is the amplitude of sinusoidal oscillations. 

The flow was perturbed by changing the direction of the ambient flow 
sinusoidally (with an amplitude of one-half of a degree) during the first cycle 
of the oscillation. The amplitude of the sine wave was the only free 
parameter. It is worth noting that this type of disturbance gradually returns 
the perturbed quantity to its initial state. 

Numerical experiments have been carried out in the range of K = 1-4, (3 
= 200, Re = 200-800, At = 0.00025-0.0005 and for various values of V r = U 0 /U m , 
as shown in Table 2. 
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Table 2 The range of the governing parameters 

(P = 200) 



V r / K = 1 2 2 4 



0.0 



X XXX 



0.60 



X 



0.65 



X 



0.70 



X 



Evidently, this is a rather limited exploration of a highly complex 
problem and requires much more numerical and experimental work. The 
purpose of the present calculations was not to provide a detailed comparison 
between the measured and calculated forces but rather to attempt to establish 
a relationship between the shedding of vortices and the relative magnitude of 
the current. 

Figures 3 through 6 show for K = 1, p = 200, V r = 0.0, and At = 0.00025, the 
in-line force and streamlines. The calculated force coefficients are Cd = 2.18 
and C m = 2.15. The corresponding theoretical values are Cd = 2.2 and C m = 2.1 
(Wang 1968). Figures 7-13 show the in-line force, the streamlines and the 
streaklines, for K = 2, P = 200, and V r = 0.0. The calculated force coefficients 
are Cd = 1-35 and C m = 2.13. The corresponding theoretical and experimental 
values are, respectively, Cd (th) = 0.95, Cd (exp) = 1.50 and C m (th) = 2.02 and C m 
(exp) = 2.0. The differences between the calculated, experimental, and 
theoretical values are due to the fact that the theoretical values do not 
account for the flow separation. Figures 14-24 show the in-line force and 
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detailed streakline plots for K = 3, P = 200, and V r = 0.0. The calculated force 
coefficients are Cd = 1.35 and C m = 2.02. The corresponding experimental 
values (Bearman et al. 1985) are Cd = 1-45 and C m = 2.10. There are no 
theoretical values to compare with for the reasons just cited. Figures 25-31 
show the in-line force and the streaklines for K = 4, p = 200, and V r = 0.0. The 
calculated force coefficients are Cd = 1.38 and C m = 1.92. The corresponding 
experimental values (Bearman et al. 1985) are Cd = 1.45 and C m = 1.95. It is 
clear from the foregoing that for the no-current case the calculated force- 
transfer coefficients are in very good agreement with those obtained 
experimentally. For the cases for which a comparison is possible, the 
calculated values also agree reasonably well with the theoretical values of 
Wang (1968). 

Figures 32-38 show the in-line and transverse forces, the streamlines and 
the streaklines for K = 4 and V r = 0.6 (with Cd= 1.17, C m = 1.65); and Figs. 39^41 
show the in-line and transverse forces and a single streakline plot for K = 4 
and V r = 0.65 (with Cd= 1.16, C m = 1-63); and finally, Figs. 42-43 show a 
streamline and a streakline plot for K = 4 and V r = 0.7 (with Cd= 1.15, C m = 
1.63). 

The calculation with current were carried out for the expressed purpose of 
substantiating the previous findings that for K = 4, within a narrow range of 
V r values, the width of the wake increases and the vortices begin to arrange 
themselves along three rows. In order to make sure that this finding was not 
a consequence of the sensitivity of the code to grid shape and size and the 
time-interval used, a more comprehensive series of calculations were 
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undertaken through the use of a non-square grid. Not only were the 
previous conclusions substantiated but also additional facts were uncovered. 
A comparison of Fig. 38 (for V r = 0.60), with Fig. 41 (for V r = 0.65), and Fig. 43 
(for V r =0.7) shows that the wake is comprised of three rows of heterostrophic 
vortices. They differ only in detail from one V r to another in the narrow 
range of V r values from about 0.6 to 0.7. At lower V r values, the inner pair of 
vortices are closer to the wake axis. As V r increases toward 0.7, the inner pair 
reaches the edges of the outer pair (Fig. 43), and the inner and the outer 
vortices form an interesting quadruple. The vortex pairs on one side of the 
street may propel themselves (as a Kelvin oval) away from the wake axis, or 
become part of a larger scale instability. In either case, it will be difficult to 
predict the behavior of such a special wake at larger times because it is rather 
difficult to distinguish between the naturally occurring instabilities and those 
occurring numerically due to truncation errors. Nevertheless, it is gratifying 
to note that the flow visualization of the early stages of the wake for K = 4, at a 
representative value of V r = 0.65, yielded results in excellent agreement with 
those presented herein. 

Table 3 shows a comparison of the calculated and experimental (Moe & 
Verley 1980) drag and inertia coefficients for K = 4 for representative values of 
V r . As expected, the inertia coefficients agree extremely well. As far as the 
drag coefficients are concerned, the agreement is not as good but certainly 
better than expected in view of the fact that [3 was 200 in the calculations and 
about 600 in the experiments. Nevertheless, the trend of the data is well 
predicted. 
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Table 3. Comparison of the Measured and Calculated Force 
Coefficients for K = 4 and Three Values of V r 



V r : 


0.60 


0.65 


0.70 


C m (ex P ) 


1.68 


1.63 


1.60 


C m <“» 


1.65 


1.63 


1.63 


Cjj (exp) 


0.93 


0.98 


1.05 


Cjj (cal) 


1.17 


1.16 


1.15 
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V. CONCLUSIONS 



The investigation reported here warranted the following conclusions: 

1. Even the higher order finite difference formulations of the governing 
equations based on the vorticity/stream-function formulation of the 
Navier Stokes equations can be solved for only relatively small 
Reynolds numbers. This is primarily due to stability and computer 
constraints. 

2. The numerical experiments with pulsating flows (oscillation plus 
steady mean flow) for K = 4 yielded force-transfer coefficients in good 
agreement with those obtained experimentally. 

3. For K = 4 and relative current velocities of about 0.6-0. 7, the vortices 
shed nearly symmetrically at each cycle and gave rise to a most 
unusual three-row vortex street, where each row is comprised of a 
pair of heterostrophic vortices. For relative current velocities larger 
than about one, the vortex wake returned to the asymmetric mode, as 
is encountered in a regular Karman vortex street. The use of a non- 
square grid in lieu of a square grid in solving the Poisson's equation 
did not alter the results. 

4. The foregoing numerical experiments could not have been possible 
had it not been due to the availability of a VAX-3520 and a CRAY 
supercomputer. It is fully realized that numerical instabilities versus 
fluid dynamical instabilities have different and at times competing 
mathematical and computational demands. It is because of this 
reason that calculations at higher Reynolds numbers (while 
preserving flow stability) and calculations at very small Keulegan- 
Carpenter numbers (while preserving numerical stability) require 
extremely large CPU times even on a supercomputer. 
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Figure 1. Grid in the Computational Domain 
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Figure 2. Grid in the Physical Domain 
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Figure 3. In-Line Force Coefficient, K = 1, Re = 200, V r = 0. 
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Figure 4. Streamlines, K = 1, R e = 200, V r = 0.0, t/T = 4.0 
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Figure 5. Streamlines, K = 1, Re = 200, V r = 0.0, t/T - 5.0 
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Figure 6. Streamlines, K = 1, Re = 200, V r = 0.0, t/T - 6.0 
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Figure 7. In-Line Force Coefficient, K = 2, Re =400, V r - 0.0 
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Figure 8. Streamlines, K = 2, Re =400, V r = 0.0, t/T = 4.0 
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Figure 9. Streamlines, K = 2, Re =400, V r = 0.0, t/T = 5.0 
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Figure 10. Streamlines, K = 2, Re =400, V r = 0.0, t/T - 6.0 
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Figure 11. Streaklines, K = 2, Re -400, V r = 0.0, t/T = 4.0 
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Figure 12. Streaklines, K = 2, Re =400, V r = 0.0, t/T = 5.0 
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Figure 13. Streaklines, K = 2, Re -400, V r - 0.0, t/T - 6.0 
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Figure 14. In-Line Force Coefficient, K = 3, Re =600, V r = 0.0 



36 




37 




38 



Figure 16. Streaklines, K = 3, Re =600, V r = 0.0, t/T = 4.0 
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Figure 18. Streaklines, K = 3, Re =600, V r = 0.0, t/T = 4.50 
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Figure 25. Inline Force Coefficient, K = 4, Re = 800, V r - 0. 
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Figure 26. Streaklines, K = 4, Re =800, V r = 0.0, t/T = 3.0 




49 




50 



Figure 28. Streaklines, K = 4, Re =800, V r = 0.0, t/T = 4.25 




51 




52 



Figure 30. Streaklines, K = 4, Re =800, V r = 0.0, t/T = 4.75 
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Figure 32. Inline Force Coefficient, K = 4, Re = 800, V r = 0.6 
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Figure 33. Transverse Force Coefficient, K = 4, Re = 800, V r = 0.6 
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Figure 34. Streamlines, K = 4, Re =800, V r = 0.6, t/T = 4.0 
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Figure 35. Streamlines, K = 4, Re =800, V r = 0.6, t/T = 6.0 
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Figure 36. Streamlines, K = 4, Re =800, V r = 0.6, t/T = 8.0 
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Figure 37. Streaklines, K = 4, Re =800, V r = 0.6, t/T = 6.0 
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TRANSVERSE FORCE COEF 
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'igure 39. 



Inline Force Coefficient, K = 4, Re = 800, V r = 0.65 




Figure 40. Transverse Force Coefficient, K = 4, Re = 800, V r = 0.65 
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Figure 42. Streamlines, K = 4, Re =800, V r = 0.7, t/T = 8.0 
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